import numpy
import sys
import random

import matrices_transpuestas as M


meses = [('ENE',(0,30)), ('FEB',(31,58)), ('MAR',(59,89)), ('ABR',(90,119)), ('MAY',(120,150)), 
               ('JUN',(151,180)), ('JUL',(181,211)), ('AGO',(212,242)), ('SET',(243,272)), ('OCT',(273,303)),
               ('NOV',(304,333)), ('DIC',(334,364))]
dictmeses = dict(meses)

Kt_MIN_MAX = [(0.031, 0.705), (0.058, 0.694), (0.051, 0.753), (0.052, 0.753), (0.028, 0.807),
            (0.053, 0.856), (0.044, 0.818), (0.085, 0.846), (0.010, 0.842), (0.319, 0.865)]


def sex2rad(grados):
    return grados * numpy.pi / 180.0

def gamma(dn):
    return 2 * numpy.pi * (dn - 1) / 365.0

def declinacion_solar(dn):
    """dn es el dia del anio, valores validos son del 1 al 365 inclusive.
    Devuelve un angulo en grados sexagecimales"""
    return (0.006918 - 0.399912 * numpy.cos(gamma(dn)) + 0.070257 * numpy.sin(gamma(dn)) - \
        0.006758 * numpy.cos(2 * gamma(dn)) + 0.000907 * numpy.sin(2 * gamma(dn)) - \
        0.002697 * numpy.cos(3 * gamma(dn)) + 0.00148 * numpy.sin(3 * gamma(dn))) * 180 / numpy.pi

#_epsilon es igual a ( distancia_media / distancia ) ^ 2, donde distancia es la distancia del sol a la tierra
# dn es un dia del anio 
def _epsilon(dn):
    return 1.00011 + 0.034221 * numpy.cos(gamma(dn)) + 0.001280 * numpy.sin(gamma(dn)) + \
        0.000719 * numpy.cos(2 * gamma(dn)) + 0.000077 * numpy.sin(2 * gamma(dn))

def irradiacion_diaria(lat, dn):
    """Calculo de la irradiacion a lo largo de un dia.
    Bo es constante solar definida como la irradiancia extraterrestre sobre una superficie normal \
    a los rayos solares y situada a una distancia del sol igual a una unidad astronomica.
    Devuelve W.h/m^2"""
    Bo = 1.367 # en Watts / m^2
    #ws es el angulo de salida del sol.
    ws = lambda lat, dn: -numpy.arccos(-numpy.tan(sex2rad(declinacion_solar(dn))) * 
                numpy.tan(sex2rad(lat)))
    #print "angulo de salida del sol: ", ws(lat, dn), "radianes"
    return (24 / numpy.pi) * Bo * _epsilon(dn) * numpy.cos(sex2rad(lat)) * \
            numpy.cos(sex2rad(declinacion_solar(dn))) * (ws(lat, dn) * numpy.cos(ws(lat, dn)) - numpy.sin(ws(lat, dn)))

def irradiacion_diaria_mensual(lat, mes):
    """Media mensual de irradiacion diaria mensual"""
    if mes not in dictmeses.keys():
        raise ValueError, '"%s" no es un mes valido' % mes
    dn1, dn2 = dictmeses[mes]
    return sum([irradiacion_diaria(lat, dn) for dn in xrange(dn1, dn2)]) / (dn2 - dn1)

# indice de claridad mensual a.k.a Ktm
def indice_claridad_mensual(lat, Gdm):
    """Gdm es Radiacion global mensual sobre superficie horizontal"""
    #Gdm = [1.99 , 2.64, 4.32, 5.32, 6.28, 7.29, 7.47,  6.62 , 5.11, 3.4, 2.16, 1.72] * 1000
    Ktm = [Gdm[i] / irradiacion_diaria_mensual(lat, mes[0]) for i, mes in enumerate(meses)]
    return Ktm

def indice_claridad(lat, Gdm):
    def hallar_matriz(valor_Ktm):
        #Este pequenio hack gracias a Arnau.
        info = [(0.3, (1, M.K_00_30)), (0.35, (2, M.K_30_35)), (0.40, (3, M.K_35_40)), (0.45, (4, M.K_40_45)),
                (0.5, (5, M.K_45_50)), (0.55, (6, M.K_50_55)), (0.60, (7, M.K_55_60)), (0.65, (8, M.K_60_65)), 
                (0.7, (9, M.K_65_70)), (1.0, (10, M.K_70_100))]
        for limit, retval in info:
            if valor_Ktm <= limit:
                return retval

    def crear_subintervalos(valor_Ktm):
        NUMERO_INTERVALOS = 10
        m = hallar_matriz(valor_Ktm)[0]
        inicio, final = Kt_MIN_MAX[m-1]
        return numpy.linspace(inicio, final, NUMERO_INTERVALOS + 1) #regresa un array
    
    Ktm = indice_claridad_mensual(lat, Gdm)
    
    Kt = [] #Lista de Kt diarios
    
    for i, mes in enumerate(meses):
        dn1, dn2 = dictmeses[mes[0]]
        for dia in xrange(dn1, dn2 + 1):
            if dia == 0:
                Kt_previo = Ktm[-1] #para 1 de enero no hay Kt_previo entonces una estima es Ktm de diciembre
            subintervalo = crear_subintervalos(Ktm[i])
            
            #determinando el subintervalo
            
            for i in xrange(len(subintervalo)): 
                try:
                    if subintervalo[i] < Kt_previo <= subintervalo[i+1]:
                        Fila = i
                except(IndexError):
                    Fila = 9   #es en realidad la fila 10 de la matriz
            MATRIZ_TRANS = hallar_matriz(Ktm[i])[1]
            
            num_aleatorio = random.uniform(0, 1)
            
            for c in xrange(11): #i son las columnas
                if MATRIZ_TRANS[Fila,:c].sum() > num_aleatorio:
                    Kt.append((subintervalo[c-1] + subintervalo[c]) / 2)
                    break
                if c == 10: 
                    Kt.append((subintervalo[9] + subintervalo[10]) / 2)
                    break
            Kt_previo = Kt[dia]
    return Kt
    
